Hard-disk equation of state: First-order liquid-hexatic transition in two dimensions 

with three simulation methods 
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We report large-scale computer simulations of the hard-disk system at high densities in the region 
of the melting transition. Our simulations reproduce the equation of state, previously obtained 
using the event-chain Monte Carlo algorithm, with a massively parallel implementation of the local 
Monte Carlo method and with event-driven molecular dynamics. We analyze the relative perfor- 
mance of these simulation methods to sample configuration space and approach equilibrium. Phase 
coexistence is visualized for individual configurations via the local orientations, and positional cor- 
relation functions are computed. Our results confirm the first-order nature of the liquid-hexatic 
phase transition in hard disks. 
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I. INTRODUCTION 

The phase behavior of hard disks is one of the oldest 
and most studied problems in statistical mechanics. It 
inspired the use of Markov- chain Monte Carlo [1] as well 
as molecular dynamics [2 j. Important progress in under- 
standing hard-disk melting [3H5] was made recently. Us- 
ing the event-chain Monte Carlo algorithm (ECMC) [6], 
a first-order liquid-hexatic transition was identified [7]. 
This transition from the low-density phase to an inter- 
mediate phase precedes a continuous hexatic-solid tran- 
sition, and thus the liquid transforms to a solid through 
an intermediate hexatic phase. 

Controversy concerning the nature of hard-disk melt- 
ing has persisted for decades. Indeed, the recently dis- 
covered first-order liquid-hexatic melting transition dif- 
fers from the standard KTHNY scenario [5HT3], which 
predicts continuous transitions both from the liquid to 
the hexatic and from the hexatic to the solid. It is also 
at variance with the first-order liquid-solid transition sce- 
nario, which exhibits no intermediate hexatic phase and 
has been much discussed [3| fT3HT8] . Near the critical den- 
sity, the system is correlated across roughly a hundred 
disk radii, and the hard-disk liquid-hexatic transition is 
thus weakly first-order [7 , with only a small discontinuity 
in density at the transition. 

For several decades, the algorithms used for this prob- 
lem [121 [I3j [T8H2Q] were unable to equilibrate systems 
sufficiently larger than the spatial correlation length to 
reliably investigate the existence and the nature of the 
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hexatic phase. This was one origin for the controversy 
surrounding this problem. Another reason was that the 
manifestation of a first-order transition in the NVT en- 
semble, and in particular the fundamental difference be- 
tween a van der Waals loop and a Mayer- Wood loop in- 
dicating equilibrium phase coexistence, was not univer- 
sally accepted in the hard-disk community, although it 
had been clearly discussed in the literature [2TH23] . 

Here, we complement and compare the recent event- 
chain results with a massively parallel implementation of 
the local Monte Carlo algorithm (MPMC) [24 , and with 
event-driven molecular dynamics (EDMD) [25]. These 
methods provide us with by far the largest independent 
data sets ever acquired for the hard-disk melting transi- 
tion, allowing us to reproduce to very high precision the 
equation of state of Ref. [7], obtain snapshots illustrat- 
ing phase separation, and compute positional correlation 
functions. 



II. SIMULATION METHODS 

A. System definition 

We consider a system of N hard disks of radius a in 
a square box of size L x L. The phase diagram of the 
system depends only on the density (packing fraction) 
r] = 7V7r<7 2 /L 2 , as the pressure is proportional to the 
temperature T. The dimensionless pressure is given by 

P* = = /3P(2<x) 2 (1) 

m(vl) 

with the inverse temperature f3, mass m and the veloc- 
ity along one axis v x . All simulations are conducted in 
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the NVT ensemble and, although our algorithms differ 
in the way they evolve the system, they all sample the 
same equilibrium probability distribution in configura- 
tion space. At finite TV, the equilibrium phase coexis- 
tence that we will observe is specific to this ensemble and 
absent, for example, in the NPT ensemble. As for any 
model with short-range interactions, the thermodynamic 
limit is independent of the ensemble. 



B. Algorithms and implementations 

Local Monte Carlo (LMC) has been a popular simula- 
tion method for hard disks [TJ [131 EH EQ] • At each time 
step, one random disk is selected and a trial move is ap- 
plied to it. The move is accepted unless it results in an 
overlap with another disk. LMC is both relatively ineffi- 
cient in sampling configuration space and inherently se- 
rial, limiting the size for which the system can be brought 
to equilibrium at high densities r] ~ 0.7 to about N ~ 10 5 
particles (see [26 for a basic discussion). Our alternative 
approaches utilize modern computer resources more ef- 
ficiently and equilibrate the system faster. They also 
provide independent checks of the equilibrium phase be- 
havior. LMC is used for comparison and as a reference 
to previous work. 

Massively parallel Monte Carlo (MPMC) [24 is a par- 
allel extension of LMC. It again applies a local trial move, 
but maximizes the number of simultaneous updates. 
MPMC extends the stripe decomposition method [27] to 
a massive number of threads using a four-color checker- 
board scheme [28]. By placing disks into cells of width 
w > 2<j, concurrent threads execute over one out of four 
subsets of cells in parallel. Within each cell, particles 
are chosen for trial moves in a shuffled order. The num- 
ber of trial moves is fixed independent of cell occupancy. 
Trial moves that would displace disks across cell bound- 
aries are rejected. The order of the four checkerboard 
sub-sweeps is also sampled as a random permutation. In 
this manner, an entire sweep over N particles satisfies 
detailed balance. The reverse sweep corresponds to an 
inverse shuffling and cell sequences with opposite trial 
moves and occurs with equal probability. To ensure er- 
godicity, the cell system is randomly shifted after each 
sweep. We implement MPMC on a graphics processing 
unit (GPU) using CUDA. Details are found in Ref. [24]. 
The MPMC simulations execute simultaneously on all 
1536 cores of a NVIDIA GeForce GTX 680. 

In event-driven molecular dynamics (EDMD), individ- 
ual simulation events correspond to collisions between 
pairs of disks [2] . The simulation is advanced sequentially 
from one collision event to the next. Between collisions, 
disks move at constant velocity (see [26] for a basic dis- 
cussion). The computation of future collisions and the 
update of the event schedule are performed efficiently 
using a binary tree and relating searching schemes [29L 
[32] . As a result, one collision event in EDMD costs only 
about ten to twenty times more CPU time than a LMC 



trial move, even for large system sizes. EDMD drives 
the system quite efficiently through configuration space 
and clearly outperforms LMC. The simulations with this 
algorithm use an Intel Xeon E5-1660 CPU with a clock 
speed of 3.30GHz. 

Event-chain Monte Carlo (ECMC) [6 replaces indi- 
vidual trial moves by a chain of collective moves that all 
translate particles in the same direction. At the begin- 
ning of each Monte Carlo move, a random starting disk 
and a move direction are selected. The starting disk is 
displaced in the chosen direction until it collides with an- 
other disk. This new disk is then displaced in the same di- 
rection until another collision occurs or until the lengths 
of all displacements add up to a total distance, an inter- 
nal parameter of the algorithm which is typically chosen 
such that the chain consists of ~ disks. With pe- 
riodic boundary conditions, ECMC is free of rejections. 
Global balance and ergodicity are preserved by moving 
in two directions only, for example to the right and up. 
ECMC is faster than LMC and EDMD [6]. The simu- 
lations with this algorithm were performed on an Intel 
Xeon E5620 CPU with a clock speed of 2.40GHz. 

For the large systems considered in this study, the ratio 
between the large absolute particle coordinates and the 
potentially small inter-particle distances becomes compa- 
rable to the accuracy of single floating point precision, so 
that round-off errors become critical. Different strategies 
allow us to cope with this problem. MPMC performs all 
computations in single-precision because today's GPUs 
run significantly slower in double-precision. We mitigate 
floating-point round-off errors by placing each particle in 
a coordinate system local to its cell. In this way, dif- 
ferences of relative positions are less affected by floating 
point precision than absolute positions. This strategy 
was also applied in [7 . EDMD calculations are per- 
formed in double-precision to fully resolve multiple co- 
incident collisions and to span the entire time domain 
from individual collision times to total simulation time. 
The ECMC algorithm is implemented for this work in 
double-precision. We use this implementation to de- 
rive high-precision numbers at the density rj = 0.698. 
Data points for ECMC at other densities are taken from 
Ref. [7 , which employs single-precision. The compari- 
son of single-precision and double-precision calculations 
at r] = 0.698 indicates that the two versions of ECMC 
yield the same result for the pressure. 



C. Pressure calculation 

The complete statistical behavior of hard disks is con- 
tained in the equation of state (pressure vs. volume or 
density), which requires the precise evaluation of the in- 
ternal pressure of the system. The equation of state al- 
lows computing the interfacial free energy and tracking 
the changes in the geometry of coexisting phase regions 
in a finite system (see [71 [2TH23] ). 

In the NVT ensemble, the pressure is a dependent ob- 
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servable. In Monte Carlo, it has to be computed from 
static configurations while, in EDMD, it may in addi- 
tion be derived from the collision rate. The disparity 
of our approaches to calculate pressure constitutes one 
more check for the implementations of our algorithms. 



1 . Pressure from static configurations 

In systems of isotropic particles with pairwise interac- 
tions, the pressure can be computed from static configu- 
rations ("snapshots") through the pair-correlation func- 
tion g(r) pp. The function g(r) is defined as the distri- 
bution of particle pairs at distance r = |x$ — Xj|, nor- 
malized such that g(r — >> oo) = 1. In practice, particle 
distances are binned into a histogram with bin size 5r. If 
n out of p pair distances are found to lie in the interval 
[r — or/2, r + or/2], then we have 



n/p 



27rrSr/V 



(2) 



The pressure P = —(dF/dV)T,N is calculated from the 
free energy F = — j3~ x \nZ and the partition function 
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where ctj = Xj/L are the particle coordinates relative 
to the simulation box. The Boltzmann weight is the 
characteristic function for overlap, i.e. zero if the config- 
uration contains overlaps and one otherwise. A change of 
volume leaves the a unchanged, but rescales the positions 
and the pair distances. 6(a\ . . . o^at) is only affected if 
one of the pair distances is at contact, hence 
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To access the contact value of g(r), we fit the histogram 
of pair distances obtained from LMC by a polynomial as 
shown in Fig. [I] and then extrapolate the fit to r = 2a 
from the right. We choose the bin size as Sr = 10 _3 a. 
The histogram is limited to r G [2a, 2.1a], and the fit is 
performed with a fourth-order polynomial. These param- 
eters are sufficient to obtain a relative systematic error 
of less than 10 -5 . In the single-precision version of the 
algorithm, g(r) shows correlated fluctuations which lead 
to systematic errors in the value of g(r) for any single bin 
(see Fig. [TJb,c)). These errors are only due to floating 
point round-off of the pair-correlation function (the sam- 
pled configurations are essentially the same). However, 
they are periodic with zero mean and do not affect the fit 
significantly. We verified that single-precision rounding 
induces a relative systematic error smaller than 10 -5 on 
the estimated value of g(2a + ). 
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FIG. 1: (a) Pair-correlation function g(r) close to contact 
for N = 512 2 at density r\ = 0.698 using LMC. Error bars 
are computed through 64 independent simulations. g(r) is 
fitted with a fourth-order polynomial. Difference between g(r) 
and the polynomial fit with the single-precision data (b) and 
double-precision data (c) for the histogram. 



2. Dynamic pressure computation 

In molecular dynamics, static configurations can be an- 
alyzed as before, but the pressure is computed directly 
and more efficiently from the collision rate via the virial 
theorem [33-35 . This avoids binning and extrapolations. 
The non-dimensional virial pressure in two dimensions is 
given by 



/3P = 



N 
V 



/3m 1 
2tZ, N 



collisions 



(5) 



where t tot is the total simulation time. The collision force 
bij = Tij • Vij is defined between the relative positions 
and relative velocities of the collision partners. In equi- 
librium, the average virial for hard disks equals [36] 



{bij) = "2a 



(6) 



Therefore, the pressure is simply given by the collision 
rate A = 1/to? the reciprocal of the mean free time to, as 



0P = 



N 
V 



1 



A 



(7) 



To test the pressure calculations, we compute the pres- 
sure with all algorithms at one representative state point. 
With each algorithm, we perform between 8 and 100 in- 
dependent runs to compute the statistical standard error. 
Results obtained with the four algorithms agree within 
numerical accuracy to < 10 -4 , which is sufficient for our 
purposes (see Table . 

III. PERFORMANCE COMPARISON 

One of the slowest processes during the time evolution 
of the hard-disk system at high density is the fluctuation 
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Algorithm #runs #disp/run f3P(2a) 2 std. error Algorithm r/#disp #disp/hour tlmc/t Speed-up 



LMC 


64 


6 x 10 11 


9.17046 


1.5 x 10" 4 


LMC 


7 x 10 5 


6.5 x 10 9 


1 


1 


EDMD 


100 


10 10 


9.17076 


1.8 x 10" 4 


EDMD 


1.8 x 10 4 


1.7 x 10 9 


39 


10 


ECMC 


32 


5 x 10 11 


9.17062 


8.7 x 10" 5 


ECMC 


2.5 x 10 4 


1.6 x 10 10 


28 


70 


MPMC 


8 


6 x 10 13 


9.17078 


4.5 x 10" 5 


MPMC 


8 x 10 5 


2.3 x 10 12 


0.9 


320 



TABLE I: Test of the pressure computations for N = 256 2 
at n = 0.698 Results of all four algorithms agree within their 
numerical accuracies. 




FIG. 2: Autocorrelation function of the global orientation 
order parameter ^(t) for N = 512 2 , n = 0.698 obtained with 
LMC, EDMD, ECMC, and MPMC. (a) Time is measured in 
number of attempted displacements (or collisions) per disk, 
(b) Time is measured in CPU or GPU hours. 

of the global orientation order parameter 

3 

which is the spatial average of the local order-parameter 
field 

1 6 

fa = -^exp(z60 j)fe ). (8) 

k=l 

The sum is over the six closest neighbors k of disk j, 
and (j)j y k is the angle between the shortest periodic vector 
equivalent to x/e — Xj and a chosen fixed reference vector. 
This approach is simpler than using the Voronoi con- 
struction [7], without affecting the autocorrelation func- 
tions. To determine the efficiency of our algorithms, we 
track the autocorrelation function of ^ [6], 

The ^6 —> ^6 + symmetry in the square box im- 
poses that C(At) decays to zero for infinite times. In 
the asymptotic limit, the decay is exponential, C(At) oc 
exp(— At/r), and we obtain the correlation time r from 
a fit of the pure exponential part (see Fig. [5| . 

We compare speeds for N = 512 2 at 77 = 0.698, that is, 
in the dense liquid close to the liquid-hexatic coexistence. 



TABLE II: Speed comparison of the four hard-disk algorithms 
for N = 512 2 , 77 = 0.698. The correlation time r is measured 
in number of displacements (disp) per disk. #disp/hour rep- 
resents the number of displacements (or collisions) per hour 
achieved in our implementations. The two rightmost columns 
show the speed-up of the algorithms in number of displaced 
disks, and in terms of CPU or GPU time in comparison to 
LMC. 

Although still a liquid, the correlation length of the local 
orient ational order parameter i/jj at this density is ~ 50<r. 
Such a large correlation length induces a long correlation 
time, and equilibration requires > 10 6 trial moves per 
disk for LMC. For the test, each algorithm is set to its 
optimal internal parameters. Total simulation run times 
are on the order of 10 3 r to 10 4 r. 

Fig. [2] illustrates that the autocorrelation functions 
decay roughly as pure exponentials. The time unit in 
Fig. [2|a) corresponds to the number of attempted dis- 
placements per disk (number of collisions in the case of 
EDMD and ECMC). As expected, MPMC decays slightly 
more slowly than LMC because trial moves across cell 
boundaries are rejected. ECMC and EDMD are sig- 
nificantly faster than LMC, confirming that these two 
methods sample configuration space more efficiently, but 
slower than MPMC. The time unit in Fig. j2^b) corre- 
sponds to the real simulation time (CPU or GPU time). 

Correlation times, number of attempted displacements 
per hour, and accelerations with respect to LMC are sum- 
marized in Table [TTJ EDMD and ECMC sample configu- 
ration space more efficiently by a factor of 39 times and 
28 times, respectively, while MPMC samples configura- 
tion space slightly less efficiently by a factor of 0.9. Ev- 
idently, the efficiency of the simple LMC algorithm is 
improved significantly with speed-ups of 10, 70, and 320 
for EDMD, ECMC, and MPMC, respectively. Although 
speeds in Table [TT] correspond to somewhat different hard- 
ware, as indicated in the methods section, the numbers 
give a clear idea of practical improvements that can be 
obtained with respect to LMC. 



IV. RESULTS 
A. Equation of state at high density 

In their seminal work, Alder and Wainwright observed 
a loop in the equation of state of hard disks [3] . As ex- 
plained by Mayer and Wood [21 (see also [22j[23]), this 
loop is a result of finite simulation sizes and therefore dif- 
fers conceptually from a classic van der Waals loop, which 
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Density rj 



FIG. 3: Equation of state from ECMC, EDMD, and MPMC 
for N = 256 2 and N = 512 2 . Error bars are mostly smaller 
than the symbols. Results agree within standard deviation. 
The inset shows the relative pressure difference AP/P of 
EDMD and MPMC with respect to ECMC for N = 256 2 . 



is derived in the thermodynamic limit. The branches of 
the Mayer- Wood loop are thermodynamically stable, but 
vanish in the limit of infinite size. 

It is known that the presence of a Mayer- Wood loop 
in the equation of state is observed in systems showing a 
first-order transition as well as systems showing a contin- 
uous transition [37] • However, the behavior of these loops 
with increasing system size is different. For a first-order 
transition, the loop is present in the coexistence region 
and is caused by the interface free energy AF. At a given 
density, the interface free energy per disk, Af = AF/N, 
can be computed by integrating the equation of state [7] . 
In two dimensions, it scales as Af oc TV -1 / 2 . In contrast, 
for a continuous transition, Af decays faster, normally 
such that AF is constant, that is Af oc TV -1 , and the 
equation of state becomes monotonic for large enough 
systems. The scaling of Af with system size, together 
with a fixed finite separation of the peaks for large system 
sizes, is a reliable indicator of the first-order character of 
a phase transition [38] . 

Fig. [3] shows the equation of state for N = 256 2 and 
N = 512 2 , obtained with ECMC, EDMD, and MDMC. 
Error bars (standard errors) are computed through in- 
dependent simulations. The inset shows the relative 
pressure difference AP/P of EDMD and MDMC with 
ECMC. Error bars correspond again to the standard er- 
ror on AP/P. We observe that the three independent 
simulations agree. In a similar way, Fig. [4] compares re- 
sults from ECMC and MPMC for N = 1024 2 . For this 
system size and the currently available computer hard- 
ware, equilibration with EDMD takes too long to be prac- 
tical. Again, all results agree within standard deviations. 
Note that while the Mayer- Wood loop shrinks with sys- 
tem size, the position of the local extrema stay fixed close 
to 77 = 0.702 and r] 0.714. 




Density 77 



FIG. 4: Equation of state from ECMC and MPMC for N = 
1024 2 . Error bars for ECMC are smaller than the symbols. 
Results agree within standard deviation. The inset shows the 
relative pressure difference AP/P of MPMC with respect to 
ECMC. 

B. Orientational order using simulation snapshots 

In [7], the orientational order parameter field was 
graphically represented for individual configurations at 
densities where coexistence occurs to show the separa- 
tion of the liquid and the hexatic phase. The character- 
istic geometry of the region of the minority phase with 
increasing density, from a circular bubble into a parallel 
stripe and again into a circular bubble, was directly ob- 
served. The spatial correspondence between variations in 
local density and orientational order provided important 
evidence for the nature of the transition. 

Ref. [7] used the projection of the local orientational 
order ipk on the global orientational order SPq and a linear 
color code. This projection is by no means a unique mea- 
sure for the orientational order parameter field. Instead, 
in Fig. [5j we use a circular color code on data obtained by 
long MPMC simulations with N = 1024 2 . At 77 = 0.700, 
the system is completely liquid and the color fluctuates 
on the scale of the correlation length. At 77 = 0.704, the 
representations once more confirm the presence of a cir- 
cular bubble of hexatic phase, visible as a large region 
of constant color (purple), whereas at r] = 0.708 a stripe 
minimizes the interfacial free energy. The hexatic phase 
is now visible in blue, while the liquid is characterized by 
fluctuations of the direction of ^q. The constant color in 
the region of the hexatic phase confirms the presence of 
the same orientational order across the system. 

C. Positional correlation functions 

To identify the nature of the high-density phase, we 
analyze the positional order in the system. The goal 
is to distinguish the hexatic phase, which has short- 
range positional order (exponential decay) and quasi- 



6 




j] = 0.700 



77 = 0.704 




0.01 



0.001 



-1/3 

77 = 0.720, ECMC - 

MPMC — 

7] = 0.718, ECMC -- 

MPMC 

°cexp(-r/120o) 



10 



100 



rlo 



1000 





270° 

77 = 0.708 ^e-phase 

FIG. 5: Snapshots of configurations obtained with the MPMC 
algorithm for N — 1024 2 . With increasing density, pure liq- 
uid (77 = 0.700), a bubble of hexatic phase (77 = 0.704), and 
a stripe regime of hexatic phase (77 = 0.708) are visible. The 
^6-phase is color-coded according to the color wheel. The in- 
terface between the liquid and the hexatic phase is extremely 
rough. 

long-range orientational order (algebraic decay), from the 
two-dimensional solid, which has quasi-long range posi- 
tional order and long-range orientational order (no com- 
plete decay). In [7], the positional order was analyzed pri- 
marily using the two-dimensional pair-correlation func- 
tion in direct space. We now complement that work by 
an analysis of the positional correlation function in re- 
ciprocal space [40] , 

Cq( r ) = (exp(zq • (x n - x m ))) . (10) 

Here, the average is done over neighboring pairs that sat- 
isfy |x n — x m | G [r — a,r + a]. With this classic method, 
the wave vector q must be chosen carefully to correspond 
to a diffraction peak. In some previous works [13 , 39 , the 
assumption was made that q would correspond to the 
reciprocal vector of a perfect triangular lattice of edge 
length ao, namely |q| = 27r/(ao y/3/2). This assumption 
is not correct because the solid phase has a finite den- 
sity of vacancies and other defects in equilibrium, which 
increases the effective lattice constant [7 . 

In the present work, q corresponds to the maximum 
value of the first diffraction peak of the structure factor 

^( q ) = ]^lZ exp ^ q ' ( Xn -x ™))- ( n ) 

n,m 



FIG. 6: Positional correlation functions C q (r) for N — 512 2 . 
Results obtained with ECMC and MPMC agree. At density 
77 = 0.718, the decay is exponential. At density 77 = 0.720, 
the decay approaches a power law oc r _1//3 . See [7] for data 
at N = 1024 2 . 



As the system can perform global rotations during the 
simulations, q rotates from one configurations to the 
other. To compute a precise average value of £(q), we 
individually rotate configurations so that is aligned 
in the same direction for all of them. 

Fig. [6] shows C q (r) for N = 512 2 at 77 = 0.718 and 
77 = 0.720. The results of ECMC and MPMC are again 
in good agreement. We observe that C q (r) decays ex- 
ponentially at 77 = 0.718. Therefore, 77 = 0.718 cannot 
be in the solid phase. At 77 = 0.720, the positional or- 
der increases drastically. C q (r) decays almost as a power 
law, r -1 / 3 , which is the stability limit for the solid phase. 
Since the coexistence phase ends at 77 ~ 0.716, the region 
77 G [0.716, 0.720] is thus hexatic and the first-order tran- 
sition observed in Fig. [3] connects a liquid and a hexatic 
phase. 



V. CONCLUSION 

We computed the thermodynamic behavior of the 
hard-disk system close to the melting transition using 
independent implementations of three different simula- 
tion algorithms to sample configuration space and two 
distinct approaches for the pressure computation. The 
equation of state data, including error bars, of Ref. [7] are 
confirmed within numerical accuracy both qualitatively 
and quantitatively. Typical relative errors are < 10 -4 , 
more than one order of magnitude smaller than finite- 
size effects for systems with up to N = 1024 2 particles. 
Such finite-size effects are manifested in the form of a 
Mayer- Wood loop in the equation of state. Our analysis 
of orientational and positional correlation functions con- 
firms the presence of a first-order phase transition from 
liquid order to hexatic order. 
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